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ABSTRACT 

We studied the statistical methods for the estimation of the luminosity function 
(LF) of galaxies. We focused on four nonparametric estimators: 1/Vmax estimator, 
maximum-likelihood estimator of Efstathiou et al. (1988), Choloniewski's estimator, and 
improved Lynden-Bell's estimator. The performance of the 1/lVax estimator has been 
recently questioned, especially for the faint-end estimation of the LF. We improved these 
estimators for the studies of the distant Universe, and examined their performances for 
various classes of functional forms by Monte Carlo simulations. We also applied these 
estimation methods to the mock 2dF redshift survey catalog prepared by Cole et al. 
(1998). We found that l/Vinax estimator yields a completely unbiased result if there is 
no inhomogeneity, but is not robust against clusters or voids. This is consistent with the 
well-known results, and we did not confirm the bias trend of l/l^ax estimator claimed 
by Willmer (1997) in the case of homogeneous sample. We also found that the other 
three maximum-likelihood type estimators are quite robust and give consistent results 
with each other. In practice we recommend Choloniewski's estimator for two reasons: 
1. it simultaneously provides the shape and normalization of the LF; 2. it is the fastest 
among these four estimators, because of the algorithmic simplicity. Then, we analyzed 
the photometric redshift data of the Hubble Deep Field prepared by Fernandez-Soto et 
al. (1999) using the above four methods. We also derived luminosity density pL at B- 
and /-band. Our B-hand estimation is roughly consistent with that of Sawicki, Lin, & 
Yee (1997), but a few times lower at 2.0 < z < 3.0. The evolution of Pl(^) is found to 
be less prominent. 

Subject headings: galaxies: evolution — galaxies: luminosity function — galaxies: 
statistics — methods: statistical 
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1. INTRODUCTION 

The luminosity function of galaxies (LF) plays a crucial role for extragalactic astronomy and 
observational cosmology. It is one of the basic descriptions of the galaxy population itself, and 
sometimes treated as a function of color (e.g. Efstathiou, Ellis, k, Peterson 1988, hereafter EEP; 
Metcalfe ct al. 1998; Lin et al. 1999) or morphology (e.g. Bingelli, Sandage, & Tammann 1988; 
Marzkc ct al. 1998), or other additional parameters of galaxies. It is also essential for interpreting 
galaxy number counts (e.g. Koo & Kron 1992; Ellis 1997) and for analyzing galaxy clustering (e.g. 
Strauss & Willick 1995; Efstathiou 1996). Furthermore, the LF is a fundamental test for the theory 
of galaxy formation (e.g. Baugh, Cole, & Prenk 1996). Recently, the exact shape of the LF has been 
of particular interest, because it is one of the key issues to the "faint blue galaxy problem" of galaxy 
number counts (Koo & Kron 1992; Ellis 1997), and may be related to dwarf galaxy formation (e.g. 
Babul k Rees 1992; Babul & Ferguson 1996; Hogg k Phinney 1997). The evolution of the LF 
is also important to derive the cosmic luminosity density, in the context of the 'Madau plot', i.e., 
cosmic star formation density as a function of redshift (e.g. Madau et al. 1996; Cowie et al. 1996; 
Sawicki, Lin, k Yee 1997, hereafter SLY97; Pascarelle, Lanzetta, k Fernandez-Soto 1998). 

Estimating galaxy luminosity function from an observational galaxy catalog is a fundamental 
work, but it is not a trivial task. Because of the flux-limited nature of the redshift survey data, 
the catalogs are inevitably censored, and suitable statistical technique is required. In the early 
stage of the extragalactic astronomy, the classical estimator, the number of galaxies in a given 
volume, $ = N/V , was used to estimate the LF (Hubble 1936). Of course this is not sufficient for 
detailed studies, and many experts have proposed ingenious methods. Schmidt (1968) invented the 
famous l/Vmax estimator in the studies of quasar population. Felten (1977) introduced the direction 
dependence of the magnitude limit. The extension for combining some different catalogs coherently 
was discussed in Avni k Bahcall (1980). Further extension to examine the evolution of the LF with 
redshift was proposed by Ealcs (1993), and the integrated variant of the Eales' estimator was used 
in the survey of the Hawaii Deep Fields by Cowie et al. (1996). Qin k Xie (1999) also developed 
this estimator with a similar line of study. The fundamental assumption of this estimator is that the 
distribution of the objects is spatially uniform. Nowadays this is regarded as a drawback, because 
we know that the galaxies have strong clustering properties in the large-scale structure. In spite of 
the drawback, 1/Knax estimator has been frequently used for extragalactic studies (e.g. Lilly et al. 
1995; Ellis et al. 1996), probably because of its simplicity in calculation. 

In order to overcome the difficulty in treating inhomogeneous galaxy distribution, some density- 
insensitive methods have been invented. Lynden-Bell (1971) proposed the C~ method, and applied 
it to the quasar sample of Schmidt (1968). This method is based on a quite sophisticated statistical 
idea as we discuss in subsequent sections. Carswell (1973) reported numerical experiments in its use. 
Jackson (1974) improved the method to combine several different catalogs, and studied the error 
estimation when the LF is expressed as an analytical form. The original method could derive only 
the shape of the probability density function, but Choloniewski (1987) (hereafter C87) improved 
the method to obtain the density normalization and to trace the density evolution simultaneously. 
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Lynden-Bell himself, and later Felten (1976) and Nicoll & Segal (1983) pointed out the drawback of 
this method that it cannot work in the faintest regime where the data points are too sparse. This 
drawback was basically overcome the introduction of smoothing method by Caditz &: Petrosian 
(1993) (CP93). Subba Rao et al. (1996) and Szokoly et al. (1998) used the method in the recent 
studies of distant galaxies. We note that this method was further generalized by Maloney & 
Petrosian (1999) to treat the doubly truncated data, but here we do not discuss it further. 

The method proposed by Turner (1979) and Kirshner, Oemler, k, Schechter (1979) used the 
ratio of the number of objects between the absolute magnitude interval [M, M + dM] and the 
number of objects brighter than M, which canceled out the density inhomogeneity. Marinoni et al. 
(1999) used this method in their analysis of the effect of the Local infall motion on the estimation 
of the LF. Similar estimators were used by Davis &: Huchra (1982) and later de Lapparent, Geller, 
&: Huchra (1989). However, as mentioned in Efstathiou (1996), this estimator does not use the 
whole sample. 

In contrast, some estimation methods using analytical LF models, which are often called 
parametric estimation methods, have been developed. Bandage, Tammann, k. Yahil (1979) (STY) 
introduced the maximum likelihood method, which was free of the effects induced by density 
inhomogeneity, in this field by using parametric Schechter form for the LF. This parametric form 
was extended for evolutionary studies of galaxies by Lin et al. (1999). Marshall et al. (1983) 
presented another parametric estimator which can treat both the LF and the evolution parameter 
simultaneously, assuming the Poisson distribution of the objects on the magnitude-redshift space. 

The maximum likelihood approach was widely used and extended to the methods which did 
not use analytical forms, often referred to as nonparametric methods. Nicoll &; Segal (1983) pro- 
posed such a type of estimator and used it for the study of their chronometric cosmology. The 
estimator which can be regarded as an advanced version of Nicoll k, Segal's method was invented 
by Choloniewski (1986) (C86). This method adopts the same assumption as Marshall et al. (1983), 
and is regarded as a binned nonparametric version of it. Another stepwise estimator, which was 
a binned analog of STY's estimator, was introduced by EEP. Now this method seems to be most 
commonly used, and is called 'the stepwise maximum likelihood method'. But note that not only 
EEP's but also most of the other estimators are based on the maximum likelihood principle. The 
EEP's method was extended to treat density evolution (Heyl et al. 1997; Springel & White 1998). 

In spite of the variety of the methods, as we see above, there had been only the comparisons of 
some methods in the literature (e.g. Felten 1976; C86; EEP; Heyl et al. 1997) before the elaborate 
intercomparison by Willmer (1997)(W97). Statistically detailed discussions are not so frequently 
seen, either, except the rigorous works of Petrosian (1992). In W97, each method was examined 
by Monte Carlo simulations and CfAl (e.g. dc Lapparent et al. 1989) data. The obtained results 
were fitted by Schechter form, and W97 discussed the distribution of the estimates by each method 
after 1000 simulations. Based on the fitting parameter distributions, W97 reported the bias trends 
for some estimators. Furthermore W97 studied the normalization estimates, and concluded that 
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the serious discrepancies between the LFs of local and distant galaxies is not attributed to the 
difference of the estimators used in the analyses. 

Now further questions arise after W97. The tests of W97 were restricted to the Schechter 
form LF. They considered, for example, the bias in the faint-end slope estimation, and concluded 
that even for the spatially homogeneous samples, estimator gives biased results. It is 

often claimed that the faint-end overestimation of the 1/Fmax estimator is caused by the density 
inhomogeneity of the Local Supercluster (e.g. Efstathiou 1996). Thus if galaxies are homogeneously 
distributed, the estimator is expected to give the correct value. If any subtler problem dwells in 
the slope estimation, further extensive experiments are required. They also mentioned the binning 
size selection. For the analysis of the recent very high-redshift data, data sparseness should be 
considered properly. 

Recently the LF of galaxies at extremely high redshift has become available with the aid of 
large telescope facilities and improved detectors. Added to this, redshift surveys have entered upon 
a new phase by development of the photometric redshift technique. The technique requires much 
lower observational cost than the spectroscopic survey, and is suitable for the analyses of the deep 
photometric data like the Hubble Deep Field (HDF; Williams et al. 1996). Though some problems 
are inherent in the technique and in the faint source finding itself (Ferguson 1998), vast advances 
have been produced by the method (e.g. Furusawa et al. 2000). The intermediate-high redshift 
results are, however, still controversial with each other. 

To settle down these problems, reliable and robust analyses of the LF are required. In this 
paper, we examined and made practical improvements for these estimation methods^. Considering 
the complicated understanding of the evolution of galaxy population, we concentrated our discus- 
sions on the nonparametric methods without any assumed functional forms for the LF. Besides we 
restricted our concerns only to the methods which use the whole sample. We used the mock catalog 
generated from various shapes of the probability density function (namely the LF). As we noted 
above, the density inhomogeneity is a basic property of the galaxy distribution. First we tested 
how accurately these methods reproduce the input density function, by using spatially homoge- 
neous mock catalogs with varying sample size. Next, we examined the estimators by using mock 
catalogs with a dense cluster and with a large void. We also used the mock 2dF catalog prepared 
by Cole et al. (1998) in this study. After checking the reliability of each method, we finally applied 
the methods to the photometric redshift catalog prepared by Fernandez-Soto et al. (1999) (FLY99) 
and studied the evolution of the LF at the very large redshift. 

This paper is organized as follows: in Section 2 we review and discuss the methods and our 

extensions. Section 3 is devoted to the tests for the performance of these methods by mock catalogs. 
We apply the methods to the photometric redshift catalog and discuss the LF evolution in section 



^Numerical calculations in this paper are based on the public software package for cosmological study 
written by one of the authors (KY). The C library can be downloaded from http://www.kusastro.kyoto- 
u.ac.jp/Tcohji/research/libcosm/. 
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4. Our summary and conclusions are presented in section 5. We briefly introduce the statistical 
model selection criterion which we used in our discussions in Appendix A. 



2. NONPARAMETRIC METHODS FOR ESTIMATING LUMINOSITY 

FUNCTION 

Before we discuss each method, we define some fundamental quantities. Let M : absolute 
magnitude, m : apparent magnitude, and di^{z) : luminosity distance in unit of [Mpc] corresponding 
to redshift z. Then 

M = m,- 5log d-L{z)- 25- K{z), (1) 

where K{z) is the i^-corrcction. Here log = log^Q. Wc use the following notation unless otherwise 
stated: 4'{M) : the luminosity function [Mpc^"^ mag^^], A/'obs ■ number of detected galaxies in the 
survey. When we use stepwise estimators, we must select the optimal binning size to suppress 
the statistical fluctuation (Sturges 1926; Beers 1992 and references therein; Sakamoto, Ishiguro, 
& Kitagawa 1986; Heyl et al. 1997). We used Akaike's information criteria (AIC: Akaike 1974) 
in order to select the optimal binning size with the least loss of information (Takeuchi 1999; for 
general discussion, see Sakamoto et al. 1986). 



2.1. Schmidt-Eales (l/Fmax) Method 

The method to construct the LF wc will discuss here was originally proposed by Schmidt (1968) 
and well known as the l/T^nax method. Bales (1993) developed it further to trace the evolution 
with redshift. Cowie et al. (1996) used this estimator in an integral form. 

We consider the absolute magnitude and redshift range 

f Ml < M < Mu 

\ zi<z<z^ ^ ' 

with a survey solid angle Q and upper and lower limiting apparent magnitude, mu and mi. Then 
we have 



J Ml 



0(M)dM = Et^-^' 

J ^ J Zmin,i 

where Zraax,i and ^min,* are the upper and lower redshift limits that a galaxy with the absolute 
magnitude Mj can be detected in the survey. We note that 

Z\ ^ -2111111,1 ■2'max,j ^ • 
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Defining z{M, m) to be the redshift that a galaxy with the absolute magnitude M is observed as 
an object with the apparent magnitude m, we get 

2max,i = min{2;u,2(Mj,mu)}, (5) 
Zrmn,i = max {zi, z(Mi, mi)} . (6) 

Actually, both Zmax,j and -Zmin.i depend on the galaxy spectral energy distributions (SEDs). Thus 
we must account for the ii'-correction when calculating the l/l^nax(i)- 

Felten (1976) proved that the Schmidt 1/T4nax estimator is unbiased, but does not yield a 
minimum variance. He also proved that the "classical estimator N/V" , which is different from the 
estimator, is biased. Willmer (1997) gave a comment that Felten (1976) had shown this 
estimator to be biased, but it is not exact. A complication of the terminology may have led to such 
a comment. 



2.2. Efstathiou-EUis-Peterson (EEP) Method 



In this subsection, wc consider the stepwise maximum likelihood method introduced by EEP. 
Since the estimator of EEP method completely cancels the density information, this method requires 
an independent estimation of the galaxy density. The EEP method uses the form of the LF 



K 



4>{M) = Y,4>kW{Mk-M) 
k=l 

The window function W(Me — M) is defined by 
W{Mi -M) = < 



AM _ AM 

1 for Mi — <M<Me + — — 



otherwise 



According to EEP ^, the likelihood function is 



iVobs 



Ef=i W{Me - Mi^ 



^{<f>k}k=i,-,Kmi} 



i/i=l,-,Afoi 



n 



1 Eti 4>iH{My,^{zi) - Me)AM 



(7) 



(8) 



(9) 



HiMn^izi) -M)={ 



1 

MliraiZi) 



AM 



+2 ^lim(^i) 







Mii^izi) - AM/2 > M 

AM/2 <M < Mii^{zi) + AM/2 
Miun{zi) + AM/2 < M 



(10) 



*Koranyi & Strauss (1997) have pointed out that the discreteness of the assumed LF causes a systematic error in 
the estimation. In order to avoid this effect, we can use the hnear extrapolated form. We do not discuss it further in 
this paper (see the appendix of Koranyi & Strauss 1997). 
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where Mii^{zi) is the absolute magnitude corresponding to the survey hmit miim at redshift Zi 
The logarithmic likelihood is expressed as 



iVn 



i=l 



J2 W{Me - Mi) In .^^ - In I E cf^eH {M^^{zi) - Me) AM \ 



(11) 



Hence, the likelihood equation becomes 



dlnC ^ ^= WjMk - Mi) _ ^= ij-(Mii,n(z,) - Mk)AM ^ ^ 
d(t>k ^ (t>k jr{Y:ti^lH{M^^{zi)-Me)AM 

and it reduces to 

A. /KM Y:7^rwiMk-M,) 



E^li ^tH{Mv,^{zi) - Mi) AM 
This equation can be solved by iteration, and we obtain the maximum likelihood estimator ^ = 

As for the normalization of the LF, some estimators have been proposed. We use the following 
estimator of the mean galaxy density, n, which was used by EEP: 



-1 -^obs -I 



where V is the maximum volume defined by the largest redshift in the sample, and ^{z) is the 
selection function, defined by 

/ (/>(M) dM 

*(z)^^%5o ■ (15) 

(/>(M) dM 



For further discussions about other normalization estimators, see the appendix of Davis & Huchra 
(1982). By combining eqs. (13) and (14), we get the final results. As Strauss &: Willick (1995) 
pointed out, at large redshift where the selection function eq. (15) is small, the estimator eq. (14) 
becomes noisy. Therefore in practice a certain cutoff should be introduced in redshift. 



2.3. Choloniewski Method 



Here we discuss the method for estimating the LF developed by C86. The advantage of the 
method is that we can obtain the density and the shape of the LF simultaneously, and can easily 
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examine the galaxy density evolution with redshift. The method explained here is an extended 
version applicable for the sample with a cosmological scale. In this subsection, we use indices i, j 
as the labels of the cells. 

We again consider the absolute magnitude and redshift range 

J Ml < M < Mu 

with a survey solid angle Q. and survey limiting magnitude miim- And let ra(r) [Mpc~^], the number 
density of galaxies in the neighborhood of the position r. If we define V as the total comoving 
volume under consideration, i.e. 



then it leads to the following expression for the mean number density as 

N 

n=y, (17) 

where N is the total number of galaxies within the redshift range zi < z < Zu- Here we adopt a 
statistical model: on the absolute magnitude-position space (M-r space) the galaxy distribution 
is /(M, r), and the probability that we find k galaxies in a volume element dMdV at (M, r), Pfe, 
is described as a Poisson distribution: 

Pu = (18) 
A = ^/(M,r)dMdF. (19) 



Here 



1 rMu 

n(r) = - / /(M,r)dM. (21) 

^ J My 



If we apply an assumption that the random variables M and r are independent, i.e. f{M,r) = 
ijj{M)u{r), then we obtain 

A = ^ ip(M) dM iy(r)dV . (22) 
n 

We integrate A over the spherical shell at redshift z and divide the M-z plane into small rectangular 
cells such that 

Mi<M <Mi+i = Mi + AM {i = !,■■■, I), ^^3) 
Zj < z < Zj+i = Zj + Az (j = 1, • • • , J). 
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Now we see that the problem is to estimate the intensity parameter Ajj inhomogeneously defined 
on the M-z plane (see Figure 1). The probability of finding kij galaxies in the cell Pkf, is 



Pk,,= , (24) 



which is characterized by the parameter 



Xij= / X = - ijjiAM lyjVj (25) 

where V[zj, -Zj+i] is the comoving volume between redshifts Zj and zj^i. Here 

I fMi+i 



I Kr)^d.. (27) 



and 



The likelihood is given by 



V, . aj ^d.. (28) 



n (29) 



and we obtain the log likelihood 



ln>C= ^ {kijluXij - Xij -hikij\} , (30) 

{Mi,Zj)&S 

where S stands for the subset of the M-z plane surrounded with Mu, Mi, Zu, 2;i, and the curve C 
defined by the selection line 

M + 5 log di^{z) + K{z) + 25 = miim . (31) 

We define the following quantities: 

w(i) = min{J,i5(j)} , jmax(^) = min{J, j5(z)} , (32) 

where Mj^Q) = {M : C n {(M,z) : z = Zj}}, and Zjg(i) = : C n {(M,z) : M = MJ} (see 
Figure 1). Using these notations, we can reduce eq. (30) as 

J Jmax (i) 

ln>C = X] X] {hj'^'^Xij - Xij -hikij^ , 

i=i j=i 

(i) 

= ^ ^ {hj In Xij - Xij -In kijl} . (33) 
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The maximum likelihood estimates (MLEs) (-00, ■ ■ ■ I'lpii^'o, - ■ ■ ji'j^'n) are the set of solutions which 
maximizes C They can be obtained, in practice, by setting the following equations to zero: 

s=l t 



g (|_iAM..v;)=o. (34, 



Sln£ /, 1 d\., d\,t 



t=l s=l 



nax ; ^ , ^ \ 



(35) 



Thus we have a set of equations which are referred to as likelihood equations. 

AM 

-^'^i 2^ z^t^t = 2^ ^it, (36) 

t=i t=\ 

^J\^ *max(i) Smax(i) 



These equations are solved by iterative procedure. At this stage, these solutions obtained here are 
not exactly the MLEs themselves, but relative values. We need one more step to obtain absolute 
values. We denote the relative solutions by '~' and exact MLEs by Prom eqs. (36) and (37), we 
have 

- V^tmax(i) I Y^tmaxCi) 7 



= " = ^^=1 ^ . (39) 

~st 



AM X^^maxO) 7 \p*max(s) , 



E^TQax 
s=l 



Then we properly normalize these solutions. Clearly it follows that 

\ij = ^M Yl Ai>jVj = Nok,. (40) 

iMi,Zj)&S iMi,Zj)eS 

If we set ipi Uj = wipi Uj, then we straightforwardly obtain the numerical factor w by eq. (40): 

w = . (41) 
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Now we obtain the LF (f){M) and density n{z) as 



1 - 



3 ' 



(42) 



AM 



n{zi) 



V 



(43) 



2.4. Lynden-Bell-Choloniewski-Caditz-Petrosian (LCCP) Method 

In this section, we discuss the method originahy introduced by Lynden-Beh (1971) as the 
'C~ method'. The estimator of this method is an analog of the Kaplan-Meier estimator used for 
censored data analyses, like survival analysis (e.g. Feigelson &; Nelson 1985; Feigelson 1992; Babu 
&; Feigelson 1996; for reference of survival analysis itself, see e.g. Kleinbaum 1996). The method 
may be the most natural application of the nonparametric statistics to the problem (e.g. Pctrosian 
1992). The rederived version of the method by C87 was improved so that it could estimate the 
LF and density evolution of galaxies at the same time. In addition, the derivation of the estimator 
was much simplified. The original method was invented to estimate the cumulative LF as a step 
function, thus the differential LF was described as a weighted sum of Dirac's (5-function. But 
obviously this form is not practical, and C87 suggested to smooth the LF. In modern statistics, 
the kernel estimator is used in the problem of nonparametric density estimation (Silverman 1986; 
Lehmann 1999). The kernel is a smooth function which is used as a substitute of the delta function, 
in order to keep the estimated density function smooth. This improvement was introduced to the 
LF estimation problem by CP93, and used for a photometric redshift catalog by Subba Rao et al. 



We unify these improvements, and show the practically convenient calculation here, which we 
call the 'LCCP method' after the names of the above contributors. We use the same notations for 

luminosity function, galaxy number density, distribution of galaxies, etc., and wc consider the same 
absolute magnitude and redshift ranges as in section 2.3. But we must note that, in this subsection, 
indices represent the labels of galaxies. This method is completely free of binning procedure. 

For the later discussion, we suppose that the galaxies are ordered as M/. < M/^^i. In the LCCP 
method, the independence assumption is also adopted for M and z, which leads to the expression 



(1996)5. 



/(M,z) = V(M)K^). 
The empirical distribution (distribution of observational data) is expressed as 



Nobs 




(44) 



k=l 



^But we note that Subba Rao et al.'s eq.(6) erroneously includes an extra exponential. 
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again N^y^s is the observed sample size, and let 

^{M) = ^^i6{M-Mi), (45) 

i=l 
JVobs 

Then the empirical distribution is 

-^obs -^obs 

/obs(M, z) = E ^» '^(^ - E '^J- - ^s{M, z) 
i=i j=i 

= E ^i^AM-Mi)5{z-Zj). (47) 
Here, xs{M-, z) is the characteristic function of the set S defined as 

Xs{M, z) = 



{M,z)iS, 

1 {M,z)eS. 



In the following discussions, the quantities M^a,^{j) and Zraa.x{i) are defined as 

-Mmax(i) = min {Mu, M5(j)} , 2;niax(i) = min {z^, zs{i)} , (48) 

where Ms(j) = {M : C n {(M, 2;) : z = Zj}}, and ^^(j) = : C n {(M, z) : M = Mi}}. Though they 
look like those used in the subsection 2.3, we note again that the indices are of galaxies. These are 
schematically described in Figure 2. Integration of eq. (47) over the interval [Mi, M^], [zj- — £, z^+e] 
(e > 0) gives 



/ / E -Me,z- ze)dzdM = 1 

J Ml Jzk-e ^^-^ 
pMu rzk+s 

= J2 ipit^jS{M - Mi)6{z - Zj)dzdM 



therefore 



A^»<^max(fe) 

1=1 



Mi<M, 



Vj E V'i = l, (i = l,---,iVobs)- (49) 



i=l 



Similarly, integration over the [M^ — £, M/j + e], [2:1, z^ gives 



V^i E ^i = l' (i = l,---,iVobs). (50) 
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Formally, we can obtain {V'j}i=i,- -,iVobs {^i}i=i>- -,^obs solving the eqs. (49) and (50), and 
the estimates of the real galaxy distribution, f{M, z), as 

/(M, -2) = ^ ^ V'i SiM - Mi) lyj Siz - zj) . (51) 

i=i j=i 



Thus 



m) = T7 f{M,z)dz=-^i^iS{M-Mi)Y,^j, (52) 
^ -^^i ^ i=i j=i 

n{z) = f{M,z)dM = — ^u,S{z-z,)J2i^^, (53) 



'A^i -• j=i i=i 

where V is the volume considered. The total number of galaxies N is 

-^obs -^obs 

i=i j=i 

These solutions are MLEs as discussed in C87. 

In spite of the clarity of the derivation, it is, actually, not an easy task to solve the eqs. (49) 
and (50) numerically if the data size N^hs is large. Thus we use the usual C estimator together, 
in order to calculate the Choloniewski's coefficients more easily. The Lynden-Bell's C~-function, 
C~(Mfe), is the number of galaxies in the region 



(55) 



I M,<M<Mk, 

\ Zl<Z< Z^^i^k) 

Let 

Cfc = C'(Mfc), fc = l,---,iVobs. (56) 
Then, using eqs. (49) and (50), we have 

i=l j=l i=l ^'^ 

and we obtain the following recursion relation: 

^fe+i = ^^r—'^k ■ (58) 
Thus, the distribution function of M (cumulative LF), $(M), is 

Ck 



$(M)cx ^ i,k=^i n (59) 



k=l k=l 



-14- 



In the real procedure, we set (Ci + 1)/Ci = 1, so the product in the above equation begins with 
k = 2. We can prove the second step of the above equation by mathematical induction. This is 
equivalent to the Lyndcn-Bcll's solution (C87). We obtain the weight {V'i}i=i,— jATobs (^^)' 
and we are able to calculate the density weight {t^j}j=i.---,Nahs ^'i- (^0). 

As we mentioned above, the weighted sum of the (5-function is not a practically useful form, and 
random fluctuation would be serious in the region where the data points are sparse. Therefore, the 
kernel estimator, which is often used in modern nonparametric density estimation, was introduced 
by CP93. This estimator is simply obtained by replacing the ^-function with a smooth kernel 
function k as 

The minimum value of the 'smoothing scale' h is restricted by the observational uncertainty, which 
was used by Subba Rao et al. (1996), but it docs not provide sufficient smoothing in general 
(CP93). The optimal value of Hm or hz may be estimated as 

hu ~ max {M^+i - Mi}-^^ ._. j^^^^ , hz ~ max {zj+i - ^j}j=i,...,jv„b, • (^1) 

It is obvious that the larger the data size N^^s is, the smaller the smoothing scale becomes. Fur- 
thermore, CP93 discussed the effect of the kernel shape on the estimates. Now it is known that 
the best shape of the kernel is parabolic, so-called the Epancchnikov kernel (Epanechnikov 1969), 
because it gives the minimum variance (Lehmann 1999; van Es 1991): 

nix) = ^il-x'). (62) 

It should be noted that, in principle, the kernel estimator is asymptotically biased, i.e. the expec- 
tation value is slightly different from the true value even if the sample size is large. 



3. TEST OF THE METHODS BY SIMULATION 

3.1. Numerical Examination with Mock Catalogs 

The validity of the estimation methods of the LF is often examined by mathematical statistics. 
For example, their statistical unbiasedness and statistical convergence were discussed in many early 
works (e.g. Felten 1976). However, quantitative evaluation frequently appears to be difficult by 
such approach, and numerical examination is quite important. Jackson (1974) used numerical 
experiments, as well as the analytical error estimation by Fisher's information matrix (see Stuart, 
Ord, &; Arnold 1999), in the study of quasar LF, and EEP also checked the errors of their method by 
Monte Carlo simulations as well as traditional information matrix approach. Mobasher, Sharpies, 
&; Ellis (1993) performed Monte Carlo error estimation to test the special method developed to 
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construct the LF at a certain waveband from the data selected at another wavelength. Heyl et 

al. (1997) examined the effect of galaxy clustering to the LF estimation by their extended EEP 
method. But computer-aided extensive intercomparison between the estimators had not been 
performed until the work of W97. They discussed the performance of several estimators when the 
LF is represented by the Schechter form 

0(M)dM = 0.41nlO</>aO-0-^("+^)(^-^*)exp (^-IQ-oAM-m,)^ ^ (g3) 

and tested the results in some cases with different Schechter parameters. Their main conclusions 
are as follows: 

1. The STY and C" methods are the best. 

2. The l/Vmax method gives biased results and tends to give higher values for the faint-end slope 
even for spatially homogeneous samples. 

3. The STY fit tends to underestimate the faint-end slope. 

4. The mean densities (normalization of the LF) recovered by most estimators are lower than 
the input values by factors (up to 20 %). 

Among these, the second one looks most strange, because as we mentioned in section 2.1, Felten has 
proved mathematically that the 1/V^nax estimator is unbiased when the homogeneous assumption 
holds. The l/l^nax method is quite frequently used in the estimation for the LFs of quasars, clusters 
of galaxies, etc., and if W97's claim is true, some widely accepted conclusions must be significantly 
changed. Thus it is necessary to examine the estimators further, not only for the Schechter form 
but for various shapes of the LF, in order to clarify the trends of the results. 

In this section, we test the four estimators discussed in the previous section by using simulated 
mock galaxy samples with a variety of the LFs which have the following functional forms: 

A. Uniform distribution, 

B. Power-law form which increases toward fainter magnitude, 

C. Power-law form which decreases toward fainter magnitude, 

D. Gaussian distribution (with standard deviation 1.67 mag), 

E. Schechter form (steep faint-end slope: a = — 1.6), 

F. Schechter form (flat faint-end slope: a = —1.1), 

with magnitude range M = [—24, —14] (Figure 3). The first three forms are designed to examine 
the effect of LF slope for estimation, and the form D, Gaussian, is to check the effect of curvature 
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of the function. Power-law LF of the form B appears ubiquitously in various types of objects. The 

form C looks apparently unrealistic, but we added this for making thorough investigation. The form 
D is interesting because approximate Gaussian form is often found in the LFs of individual galaxy 
types. We applied Box-Muller method (Box & Muller 1958) to generate Gaussian distribution 
from uniform random number, and von Neumann's acception-rejection method to obtain other 
distributions (see Knuth 1998 for details). We set the sample sizes ~ 100 and ~ 1000, to study the 
behavior of the statistical estimators with galaxy number. Here 'sample size' means the detected 
number of galaxies after magnitude selection (observation) procedures. Therefore the underlying 
population density for each LF form is different from each other. The estimation of galaxy spatial 
density is an important part of the derivation of the LF. What to be estimated is the total galaxy 
number including the galaxies too faint to be observed. In our simulations, we stochastically 
produced galaxies according to the assumed LF, distribute them in space, calculate their observed 
flux, and judge that they could be observed or not. Therefore, the total number corresponds to the 
number of Monte Carlo trials. We fixed the number of trials through one sequence of simulations 
with a certain LF shape and spatial density. 



3.1.1. Mock Catalog with Spatially Homogeneous Distribution 

First we construct a set of mock galaxy samples with spatially homogeneous distribution in 
order to investigate the bias trend of the estimators, especially for the Schmidt and Eales' 1/Fmax- 
We set the redshift range up to 0.1, and we adopted the Hubble parameter Hq = 75 kms~^Mpc~^, 
Oo = 0-2 {qq = 0.1), Ao = 0, and limiting magnitude miim = 13 mag in the series of simulations. No 
i^-correction is considered here. We constructed 100 representations for each LF form and sample 
size, and applied the four estimators to each sample. 



3.1.2. Mock Catalog with a Dense Cluster and with a Void 

We, next, investigate the response against density inhomogeneity of galaxies. We consider some 
extreme cases for clear understanding. For the case with density enhancement, we constructed a 
series of mock catalogs with a dense spherical clump, to which half of the galaxies belong. The 
clump lies at a distance of 0.8 Mpc, and its radius is 0.8 Mpc. We call the mock catalog the "cluster 
sample" . An example of the spatial configuration of galaxies of a cluster sample is described in 
Figure 4. Then we also constructed a set of the mock catalogs with a large spherical void without 
galaxies. The void lies at a distance of 0.8 Mpc and its radius is 1.6 Mpc. We call this mock 
catalog the "void sample" . The overall underlying density of cluster and void samples defined in 
a considered volume is the same as the homogeneous samples for each LF shape, i.e. we set the 
number of Monte Carlo trials the same as that of the homogeneous sample for each LF shape. 
Therefore, the observed sample size of the cluster sample is larger than that of the homogeneous 
sample, because we put the dense clump in the considered volume. In the case of the void sample. 
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the observed galaxy number is smaller than that of the homogeneous one. 

3.1.3. Results 

The results for the 1000-galaxy samples are shown in Figure 5 - 10. The solid lines represent 
the input distributions, and the symbols are the averages of the estimates. The error bars depict 
the standard deviations of the mean of the estimates for 100 representations. Figures 5a, 5b, and 
5c are the results from the spatially homogeneous sample, from the cluster sample, and from the 
void sample, respectively. This is also the same for Figures 6 - 10. 

At a glance, we see that all estimators give consistent results with each other, and we do not 
find any bias trends in our numerical experiments for any LF forms in the case of homogeneous 
samples. For cluster samples, the 1/V^ax method yields strongly distorted estimations, as widely 
recognized. The overestimation corresponding to the dense clump clearly appeared in the 1/P^ax 
results. In contrast, the other three estimators were not affected by the dense cluster at all. The 
estimates appeared to be consistent with each other, and showed perfect agreement with the input 
LFs. The 1/T4nax method was also affected by the large void, and gave underestimated results. 

Large fluctuations appear at the faint end of the LF, because the number of available data 
points is small, especially in the case of the LF form C and D. We can obtain statistically stable 
estimates if the slope is properly steep, and the more shallow the slope is, the larger the fluctuation 
becomes. This is clearly shown in Figure 5 - 10. 

In principle, the error bar of the Choloniewski method is larger than those of the other methods, 
because the method subdivides the M — z plane both in M and z. This procedure enables us to 
estimate the shape, the normalization, and the evolution of the LF at the same time. On the other 
hand, this becomes a drawback when the data size is small, because the shot noise dominates. 
Therefore we cannot expect a firm estimation with the Choloniewski method when the sample size 
is smaller than 100. 

Here we mention the calculation time that each method needed for the same sample size. 
Because of its algorithmic simplicity, Choloniewski method is the fastest among the four methods. 
When we analyze the 1000-mock data, the relative calculation times of the l/V^ax; EEP, and LCCP 
methods normalized with that of Choloniewski method are 2.76, 2.73, and 1.87, respectively. This 
advantage is quite significant when we treat a large sample of ~ 10^~^ galaxies. We estimate the 
LFs from large datasets of sample size 250,000 in Section 3.2 by the l/Vmaxj EEP, and Choloniewski 
methods. The relative calculation times of the 1/V^ax and EEP methods normalized with that of 
Choloniewski method are, in this case, 15.01 and 131.74, respectively. The EEP method takes 
longer calculation time because it needs more iterations in the procedure than others do. The 
1 / Knax method derives the maximum volume Knax for each galaxy, and also needs some calculation 
time. The LCCP method requires a large stack for data sorting procedure, which is a requirement 
of this method. Thus we stress that the Choloniewski method is most economic from the standpoint 
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of practical computing. 

Figures 11-16 are the same as Figures 5 - 10, except that the data size is 100. We see it 
is often not possible to determine the faint end of the LF accurately for such small datasets. The 
fluctuation became larger than the result of the 1000-samplc, but we did not find the systematic 
bias trend from our results. Thus we conclude that when the galaxy distribution is homogeneous, 
all four estimators provide the consistent and correct results, even the 1/V^max estimator. 

3.2. Mock 2dF Redshift Catalog 

The Anglo-Australian 2-degree field (2dF) galaxy redshift survey is now underway^. This 
survey will measure 250,000 redshifts, up to 2 ~ 0.2, and be complete to an extinction corrected 
apparent magnitude of 6j < 19.45 mag. In order to develop statistical methods and faster algo- 
rithms for the analyses of such large upcoming redshift surveys. Cole et al. (1998) prepared an 
extensive set of mock 2dF catalogs constructed from a series of large cosmological A/"-body sim- 
ulations. The simulations span a wide range of cosmological models, with various values of the 
density parameters, Qq, the cosmological constant, Aq, and the shape parameter T and amplitude 
of the density fluctuation erg. The LF is assumed to be a Schechter form with the parameters 
reported by APM-Stromlo bright galaxy survey (Loveday et al. 1992), M^j* — 51og/i = —19.5 
mag, a = —0.97, and ^* = 1.4 x 10~^ Mpc""^. The K-correction is assumed to be canceled by 
evolutionary correction. 

We applied the three methods to the mock 2dF catalog in order to see how accurately they 
can reproduce the true LF when they are used in the analysis of realistic large redshift surveys. We 
did not use the LCCP method for this sample. When we treat such a large catalog, the advantage 
of the Choloniewski methods is extremely significant. We also focused on the difference between 
the real-space data and the redshift-space data which is affected by the redshift distortion. The 
redshift distortion causes a scatter in the estimated luminosities of galaxies. In this study, we used 
three mock catalogs, named El (Einstein-de Sitter: = l,Ao = 0,r = 0.5, as = 0.55), L3S 
{Qo = 0.3, Ao = 0.7, r = 0.25,0-8 = 1.13), and 038 (Qq = 0.3, Aq = 0,r = 0.25,^8 = 1.13). The 
catalogs we selected are all cluster- normalized, i.e. the amplitude of the initial power spectrum 
is set to reproduce present abundance of rich galaxy clusters in the local Universe (e.g. Viana & 
Liddle 1996; Kitayama & Suto 1997) and /i = F/l^o- 

We compare the input LF and the estimated LF in Figures 17, 18, and 19. Figure 17 shows 
the LF derived from the Einstein-de Sitter (EdS) Universe, Figure 18 is the LF derived from L3S 
data, and Figure 19 is the LF derived from 03S data. The left panels in these Figures show the LFs 
derived from the redshift-space data, and the right panels, those from the real-space data. First, 
we see that all the estimators provided perfectly consistent results, and they show an excellent 



See http://msowww.anu. edu.au/'colless/2dF for recent status. 
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agreement with the input LF. There are no significant difference between the real- and redshift- 
space datasets. The sUght deviations of l/I^nax estimates are caused by the clustering in the 2dF 
mock catalog. Thus wc do not have to consider the redshift distortion effect seriously when we 
derive the galaxy LF from such large-volume redshift surveys. When we use such a large survey, 
we should rather mention the photometric calibration as a more important error source. 

4. APPLICATION TO THE HUBBLE DEEP FIELD 

Recently some authors claim that the faint-end slope of the LF becomes steeper with redshift 
at z < 1 (e.g. Ellis et al. 1996; Heyl et al. 1997; but see Lin et al. 1999). The LFs for some special 
classes of galaxies such as Lyman-brcak objects (Steidel et al. 1998) or Ly-a emitters (Pascarelle 
et al. 1999) are now also available. Wc, however, do not have a coherent understanding of the 
evolution of the LF and the evolution of the luminosity density, pL- At low redshift, Zucca et 
al. (1997) reported a high normalization LF with cj)^ = 0.020 /i^ Mpc~^, and Ellis et al. (1996) 
obtained cf)^ = 0.026 /i^ Mpc~^, while Loveday et al. (1992) derived = 0.014 /i^ Mpc"^, and Las 
Campanas Redshift Survey result (Lin et al. 1996) is similar to the value of Loveday et al. (1992). 
The local value of the LF parameters plays a crucial role in the study of galaxy evolution, since it 
controls the redshift dependence of pL- Cowie et al. (1999) showed a rather mild evolution of the 
UV luminosity density at z < 1 from their surveys. On the other hand, high redshift LF estimations 
are also controversial with each other. Gwyn & Hartwick (1996) claimed dramatic changes in the 
LF from z = to z ^ 5, becoming fiat between — 24 < Mb <— 15for3<2;<5. On the contrary, 
SLY97 reported more familiar Schechter form with a = —1.3 for the LF at 3 < z < 4. Mobasher 
et al. (1996) suggested a stronger evolution of the LF. From a deep multiband photometric survey, 
Bershady et al. (1997) gave a constraint which ruled out Gwyn & Hartwick's result. 

Thus in this section, we apply the four estimators to the photometric redshift catalog of the 
HDF to study the evolution of the LF shape. For the observational data, the error estimation is 
complicated, because the estimation procedure of the LF involves the magnitude selection, weight- 
ing, etc. In such cases, bootstrap resampling analysis is known to be often superior to classical 
analytic methods in order to estimate statistical properties (e.g. Efron &; Tibshirani 1993; Babu & 
Feigelson 1996; Davison &; Hinkley 1997). Thus we used the bootstrap method for the estimation 
of the statistical uncertainties. When we perform the bootstrapping, how to generate good random 
numbers is important. We generated the uniform random number by Mersenne Twister method^ 
developed by Matsumoto & Nishimura (1998). 



^For recent development, see http://www.math.keio.ac.jp/matumoto/mt.html. 
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4.1. Sample 

We used the photometry and photometric redshift catalog of the HDF prepared by FLY99. 
Their catalog contains 1067 galaxies, with AB(8140) < 26.0. The photometric redshifts are derived 
based on both UBVI (F300W, F450W, F606W, and F814W, respectively; Williams et al. 1996) 
obtained by WFPC2, and Ji/if obtained by the IRIM camera on the Kitt Peak National Observa- 
tory 4-m telescope. The object detection and photometry are performed using SExtractor (Bertin 
& Arnouts 1996). Details of the procedures are found in FLY99. In the peripheral region of the 
WFPC2 image (referred to as zone 2), the detection limit is AB(8140) = 26 mag, and in the inner 
region (zone 1), AB(8140) = 28 mag. We restricted our analysis to the inner zone 1 sample. The 
solid angle of zone 1 is 3.92 arcmin^ = 3.32 x 10~^ sr. The sample size is then 946 galaxies. 

FLY99 used four spectral templates given by Coleman, Wu, k. Weedman (1980) to determine 
the photometric redshifts. For ultraviolet wavelengths, the templates are extrapolated by using 
the results of Kinney et al. (1993), and for infrared, by the models of Bruzual k. Chariot (1993). 

Evolutionary corrections are not included in the model spectra to avoid additional parameter depen- 
dence. According to Coleman et al. (1980), they classified the galaxy spectra into four categories: 
1. Elliptical, 2. Sbc, 3. Scd, and 4. Irr. We used these labels to set the i^-corrections. 

In principle, the SED must be the same as the templates used in FLY99, but for simplicity and 
comparison with other studies, we used the galaxy SED sample compiled by Kinney et al. (1996). 
The data of Kinney et al. (1996) have almost the same properties as the SED templates of FLY99, 
thus we can use them comfortably. To construct the K-correction function, we first selected the 
sample galaxy SEDs corresponding to the labels of FLY99, and fitted polynomial functions from 
1st order to 6th order. The order of polynomial fitting was decided by referring to AIC, and we 
chose the 5th order. 

4.2. Results and Discussions 

We show the redshift-dcpcndent LF at /-band and B-hanA in Figures 20 and 21, respectively. 
The symbols represent the estimated LFs by the four methods. We show the LFs of the HDF 
at < z < 0.5 (106 galaxies), 0.5 < z < 1.0 (193), 1.0 < z < 1.5 (204), 1.5 < ^ < 2.0 (193), 
2.0 < z < 3.0 (117), and 3.0 < z < 6.0 (109). The sample is /-band selected, and we derived the 
5-band LF by following the discussion of Lilly ct al. (1995). We stress that the four different LF 
estimators give consistent results for the HDF sample, same as the results for the mock catalogs. 

We clearly see the evolutionary trend of the LF with redshift. But we note that, though we can 
fit Schechter function, it is not so easy to derive the parameters a, M*, or 0* precisely, because the 
Schechter function is rather smooth and the errors of these characteristic parameters are strongly 
correlated. These parameters can be easily affected by statistical fluctuations. We will discuss more 
details of the /- and B-band results at each redshift range in the following. 
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4-2.1. I-band LF at < z < 0.5 

In Figure 20, the dotted line represents the local I-band LF obtained by Metcalfe et al. (1998). 
Metcalfe et al. (1998) pointed out a possible upturn of the faint end of their multiband LFs, though 
they took a prudent attitude in concluding firmly. The upturn magnitude Mj ~ — 15 + log/i mag 
(in Figure 20, h = 0.75) is in good agreement with that of our lowest redshift LF except the 
normalization. 

4.2.2. Evolution at B-band: < z < 0.5 

We compare the normalization of the LF with other previous results. Our B-hand LF shows 
roughly good agreement with other local LFs. In Figure 22 wc put our LF, SLY97 Schechter fit, a 
nd Schechter functions reported by Metcalfe et al. (1998) and Ellis et al. (1996) (Autofib Redshift 
Survey) The dotted line depicts SLY97 LF, dot-dashed line represents the Metcalfe et al. (1998) 
B-band LF, and long dashed line is the Autofib LF at z < 0.1. Our LF and that of SLY97 agree 
with higher-normalization LF reported by Autofib Survey, but are significantly higher than that of 
Metcalfe et al. (1998) Autofib LF is also consistent with the LF of ESP Survey (Zucca et al. 1997), 
while Metcalfe ct al.'s LF is consistent with EEP LF and Stromlo-APM LF (Lovcday ct al. 1992). 
But since the solid angle covered by HDF is extremely small and thus the normalization can be 
strongly affected by cosmic variance, we should not go into further discussion. 

We note that the SLY97 M* value is significantly higher than those of other surveys. This is 
because the exponential decline at bright end is not observed in the HDF LF at < z < 0.5, and 
a bump exists at Mb ~ — 20 mag. We should also mention that the error bar of the of SLY97 
is very large (1.6 mag). Considering the large error bars and the uncertainty of the photometric 
redshift, we conclude that the bright M* is not a real feature. 

At this lowest redshift, the rise of the faint end is prominent. The problem of the faint-end slope 
of galaxy LF has long been a matter of debate, and we do not have a widely accepted consensus 
yet. As we already pointed out in the above, even in /-band we find a steepening of the faint end. If 
this steep faint end is the artifact of the clustering, the LF derived from 1/Knax and those derived 
from other estimators should have been different (Section 3.1). But in fact, they are consistent 
with each other. Thus we conclude that, at least in the HDF, the faintest end of the LF has a steep 
slope in the Local Universe. 

4.2.3. Evolution at B-band: 0.5 < 2: < 1.0 

It seems that the brighter galaxies are more numerous than the local value at 0.5 < z < 1.0. 
Here we should remember the fact that the "fuzzy" redshift determination is known to affect the 
shape estimation (Liu et al. 1998). Liu et al. (1998) showed by numerical experiments that the 
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faint-end slope is underestimated and is overestimated by the photometric redshift blurring. 
The uncertainty of the photometric redshift is rather independent of the object redshift, so the 
effect will be severer at the low-z, and the can be overestimated. Thus the increase of the 
bright galaxies is partially due to this effect. But we can discuss the trend of the LF evolution by 
comparison of the LF derived from photometric redshifts consistently (Liu et al. 1998). 

4.2.4. Evolution at B-band: 1.0 < z < 2.0 

The LFs of the redshift range 1.0 < 2; < 1.5 and 1.5 < 2; < 2.0 are the most reliable ones among 
the LFs in Figure 21, since the sample size is twice larger than those of the other redshift ranges, 
and in addition, the photometric redshift error becomes worse again at 2; > 2. SLY97 suggested 
the steepening of the faint-end slope at this redshift. Our LF of 1.5 < z < 2.0 presents a similar 
feature, though the slope becomes flatter at the faintest regime. The deformation of the LF from 
2 ~ to 2 ~ 2 supports that the steepening of the fainter side of the LF, which is confirmed at 
2; < 1, is continued up to 2; ~ 2. We do not find a significant shift of at this redshift range. 

4-2.5. Evolution at B-band: 2.0 < z 

The normalization of the furthest redshift LFs settles down to the local value, while we also 
find a brightening of M* at 2; > 3.0. We must be careful that in such high redshift, cosmological 
surface brightness dimming is quite severe, and selection effect becomes significant (Ferguson 1998; 
Weedman et al. 1998). Other kinds of selection effects are discussed in Pascarelle et al. (1998). 
Therefore, there can exist more numerous galaxies than estimated. Further discussions require 
delicate treatment of such effects. 



4.2.6. Luminosity density evolution 

In order to explore the cosmic star formation history, we derived the luminosity density at B- 
and /-band based on our LFs. We fit the Schechter function and extrapolate the faint end below the 
detection limit. As we mentioned above, the Schechter parameters are poor indicators of the galaxy 
evolution, but the integrated luminosity density pl is regarded as an indicator of the evolution of 
galaxies, because in case the Schechter parameters are significantly affected by the fluctuations, pL 
is robust against the effect. We showed the derived pL in Figure 23. The upper panel shows the 
evolution of the S-band luminosity density, pi,{B), and the lower, the /-band luminosity density, 
Pl{I)- Open squares are Pl{B) derived from CFRS (Lilly et al. 1995), open circles, pi,{B) from 
Autofib (Ellis et al. 1996), open triangle represents the value from Stromlo-APM (Loveday et al. 
1992) and open diamond, ESP value (Zucca et al. 1997). Crosses are the estimates of SLY97. In 
this paper we did not try to correct for the reddening effect of dust. 
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We see the local diversity of the Pl{B), corresponding to the normalization discrepancy in 
Figure 23. Despite the fact that the local LF is hard to derive from the HDF data, our \ow-z value 
is consistent with other previous results. Added to this, our Ph{B) at 0.5 < z < 1.0 significantly 
suffers from the redshift blurring effect, but it is also consistent with CFRS highest redshift point 
within the errors. As a whole, our result is consistent with that of SLY97, except for 2.0 < ^; < 3.0. 
In this redshift range, Pl(-B) of SLY97 is several times larger than our estimate. This difference 
may be because SLY97 obtained a steeper a and brighter than ours. We find a flatter LF 
slope, and the estimates fainter than —20 mag arc not reliable in our result since the fluctuation is 
horribly large at this redshift. If we choose steeper slope, our pL will be higher. We have to wait for 
larger datasets to address this problem. At very high-z, we derived a moderately high pL, implying 
significant numbers of stars have already formed at such a high redshift. The evolution of pl(^) 
appears to be flat. At the longer wavelength, the observed light is dominated by the contribution 
from lower-mass stars, and the temporal change of the SFR is less prominent. We, in addition, 
should note that the /-band results are subject to larger i^-correction extrapolation uncertainties, 
compared with the 5-band results. 

At last, we must notice that the above discussions do not account for the fact that the sample 
is selected at /-band, and the selection criterion is different for each redshift range. At ^ < 0.5, the 
sample is safely regarded as /-selected, while at z > 3.0, they are in fact rest UV-selected. Thus the 
ideal discussion on the evolution of the galaxy LF should be based on the suitably designed survey 
as performed by Cowie et al. (1999). We will consider this point, and make more sophisticated 
discussions elsewhere (Takeuchi 2000, in preparation). 

5. SUMMARY AND CONCLUSION 

The estimation of the LF from observational data is not a trivial task, because of the flux- 
limited nature of the astronomical data. We focused on the following four estimators: 1) Schmidt- 
Eales (1/Kiiax) method , 2) Efstathiou-Ellis-Peterson (EEP) method, 3) Choloniewski method, 
and 4) Lynden-Bell-Choloniewski-Caditz-Petrosian (LCCP) method. We improved some of the 
estimators for studying the very distant universe, and examined their performances for much wider 
class of functional forms by Monte Carlo simulation. We tested these four estimators by the 
numerical experiments with mock catalogs. We also used the mock 2dF catalogs prepared by 
Cole et al. (1998). Then we applied these estimators to the HDF photometric redshift catalog of 
Fernandez-Soto et al. (1999). Our conclusions are as follows: 

1. If the sample is spatially homogeneous, all estimators give consistent results with each other, 
and we did not find any bias for any LF shapes. Thus, when we have a sufficiently large galaxy 
sample, we can use any of the estimators examined in this paper. Even when the sample size 
is smaller, the mean values remain unbiased, though the standard deviations become larger. 

2. Large fluctuation appears at the faint end of the LF, because the amount of available data is 



-24- 



small. Therefore, the flatter the LF slope is, the larger the fluctuations become. When the 

sample size is small, fluctuations in the Choloniewski method become seriously large due to 
shot noise, and thus we recommend this method for the analysis of large samples. 

3. When a large cluster or void exists, l/14iax estimator is severely affected in its LF shape 
estimation. The other three estimators are not affected by a cluster or void at all. They 
gave consistent results with each other, and the estimates showed perfect agreement with the 
input LFs. 

4. We examined the calculation time of each method. Because of its algorithmic simplicity, 
Choloniewski method is the fastest among the four methods. The EEP method needs more 
iterations in the procedure than others do, and longer calculation time. The l/Knax method 
calculates the maximum volume V^ax for each galaxy, and also needs significant calculation 
time. The LCCP method requires a large stack for data sorting procedure, which is a require- 
ment of this method. Thus we stress that the Choloniewski method is the most economic 
from the standpoint of practical computing. 

5. We examined more realistic large mock samples, specifically mock 2dF catalogs prepared by 
Cole et al. (1998). We found that the redshift distortion does not affect the LF estimates. 
When we treat such a large catalog, the advantage of the Choloniewski method is extremely 
significant in terms of the computation time. 

6. We derived the I- and B-band luminosity function of the HDF. The four different LF esti- 
mators gave consistent results for the HDF sample. We found the overall brightening of the 
LF. It seems that the faint-end steepens toward z = 2 — 3, and settles down to the local 
value at 2; ~ 3. We note that the "fuzzy" redshift determination is known to affect the shape 
estimation (Liu et al. 1998). 

7. We found a rather mild evolution of the LF. Despite the fact that the local LF is hard to 
derive from the HDF data, our low-2; value is consistent with other previous results. Our 
Pl{B) at 0.5 < 2; < 1.0 is also consistent with CFRS highest redshift point within the errors. 

As a whole, our result is roughly consistent with that of SLY97, but lower at 2.0 < z < 3.0. 
At very high-z, we derived a moderately high pL, implying that a significant numbers of stars 
have already formed at such a high redshift. We found that the evolution of Pl{I) is flat. 
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A. AKAIKE'S INFORMATION CRITERION 

In this appendix, we make an informal introduction of Akaike's theory. The meaning of the 
maximum likelihood method is clearly understood by using the concepts of information theory. 
Since the middle of 1970's, vast advances have been made in the field of the statistical inference by 
the discovery of Akaike's Information Criterion (AIC: Akaike 1974). The AIC is closely related to 
the information entropy, especially to the 'relative entropy' of two probability distributions. The 
relative entropy has a property just like a distance in differential geometry, i.e. it is a distance 
between the two probability distributions. Using AIC enables us to compare the goodness of a 
certain model with that of another type directly. For this fascinating property, AIC is applied to 
various fields of studies. The AIC is expressed as 

AIC = -2{\nC{9)-K) , 

where C is a likelihood function, is a set of maximum likelihood estimators, and K is the number 
of free parameters of the assumed model. The "most preferred" model is the one which minimizes 
the AIC. 

Here we present the problem of polynomial regression model selection by using AIC as an 
example. As we mentioned in Section 4, we adopted this procedure to determine the order of 
isT-correction as a function of redshift. Given a set of n pairs of observations {xi,yi), • • ■ , 
we fit the m-th order polynomial model 

m 

where is an independent random variable which follows the normal distribution with mean 
and dispersion o"^. The variables Xi and y-i are called the explanatory variable and the objective 
variable, respectively. This model is a conditional distribution of which the distribution of the 
objective variable y is a normal distribution f{yi) with the mean oq + aiXi + • • • + a-mX^ and the 
variance cr^, i.e. 



/(yi|ao, • • • , a^, g^) = ^=== exp < 



^(y,-Y.atx\\ 

\ ^=0 / 



(A2) 

Therefore, when a set of data is (a;i,?/i), • • • , {xn,yn), the likelihood is given by 

n 

J^iUi,- ■ ■ ,yn\ao,- ■ ■ ,am,o-'^) = JJ/(yi|oo,---,am,£^^^ 

i=l 
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1/ \ l 



The log likelihood is then expressed as 

^ n / m \ 2 

ln/:(y|ao,---,am,f^^) = - 2 ^^^^^ ~ 2 ^^^^ ~ 2^ ^ ( ~ ^ '^^^N ' ^^^^ 

''^ i=i \ e=o ) 

The log likelihood eq. (A4) is maximized with respect to oq, • ' ' > o^m when 

n / m \ ^ 

i=l V £=0 / 

is minimized. Thus, in the case of polynomial model fitting, the maximum likelihood procedure is 
equivalent to the least square method. The necessary conditions that ao, • • • maximize S are 
the normal equations of the least square. 



dS 
dao 

dS 

dai 

dS 
dam 



2^ =0 

1=1 \ £=0 J 
n / m \ 



-2j]xr 2/.-^a,xn = 0, (A6) 
1=1 \ e=Q J 



and the maximum likelihood estimates oq, • • • are obtained by solving these linear equations. 
Besides, the necessary condition that maximizes eq. (A4) is 



5 ln£ n 



d 2(72 

The maximum likelihood estimate of the residual variance is 

i=l \ ^=0 / \i=l i=0 1=1 / 

Hereafter we denote the residual variance for a model with m-th order as a^{m). Then, from 
eqs. (A4) and (A8), the maximum log likelihood becomes 

ln>C(y|ao,- • ■ ,am,(T'^) = --ln27r- -lncr2(m) - - . (A9) 
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The m-th order polynomial model has m + 2 parameters (oq, • • • ,0^,(7^(771)). Substituting K = 
m + 2 and eq. (A9) into eq. (Al) gives the AIC of the m-th order model, 

AlC(m) = n(ln27r + l)+nlna2(m) + 2(m + 2) . (AlO) 

The result of our polynomial fitting to the if-correction is presented in Figure 24. We also 
summarize the AIC value for each polynomial order in Table 1. The AIC values of elliptical, Sbc, 

and Irr in table 1 really took their minima in the case that the fitting polynomials were those with 
5th order, and only Scd data preferred the 6th order. Putting all accounts together, we chose 5th 
order polynomial model. 
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Figure Captions^ 

Fig. 1. — Schematic description of the Choloniewski's method. 

Fig. 2. — Schematic description of the Lynden-Bell-Choloniewski-Caditz-Petrosian method. 

Fig. 3. — Simulated functional forms for the luminosity functions. The density scale is arbitrary 
because it depends on the adopted sample size, and we set the scale for the case of size ~ 1000. 

Fig. 4. — An example of the spatial configuration of galaxies of a cluster sample. 

Fig. 5a. — Estimates of the luminosity function (LF) for mock samples with size ~ 1000. The input 
luminosity function is with the functional form A. This is the LF of the spatially homogeneous 
sample. 



Fig, 


. 5b.- 


- The 


same 


as 


Figure 


5a, except that this is the LF estimated from the cluster sample. 


Fig, 


. 5c.— 


- The 


same 


as 


Figure 


5a, except that this is the LF estimated from the void sample. 


Fig, 


. 6a. — 


- The 


same 


as 


Figure 


5a except that the LF shape is B. 


Fig, 


. 6b.- 


- The 


same 


as 


Figure 


5b except that the LF shape is B. 


Fig, 


. 6c. — 


- The 


same 


as 


Figure 


5c except that the LF shape is B. 


Fig, 


. 7a.— 


- The 


same 


as 


Figure 


5a except that the LF shape is C. 


Fig, 


. 7b.- 


- The 


same 


as 


Figure 


5b except that the LF shape is C. 


Fig, 


. 7c.— 


The 


same 


as 


Figure 


5c except that the LF shape is C. 


Fig, 


. 8a.— 


- The 


same 


as 


Figure 


5a except that the LF shape is D. 


Fig, 


. 8b.- 


-The 


same 


as 


Figure 


5b except that the LF shape is D. 


Fig, 


. 8c.— 


The 


same 


as 


Figure 


5c except that the LF shape is D. 


Fig, 


. 9a.— 


- The 


same 


as 


Figure 


5a except that the LF shape is E. 


Fig, 


. 9b.- 


' The 


same 


as 


Figure 


5b except that the LF shape is E. 


Fig, 


. 9c.— 


The 


same 


as 


Figure 


5c except that the LF shape is E. 



Fig. 10a. — The same as Figure 5a except that the LF shape is F 



^AU figures listed below are available from ftp:/ /ftp. kusastro.kyoto-u.ac.jp/pub/kohji/lf/ . 
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Fig. 10b. — The same as Figure 5b except that the LF shape is F. 
Fig. 10c. — The same as Figure 5c except that the LF shape is F. 

Fig. 11a. — Estimates of the luminosity function (LF) for mock samples with size 100. The input 
luminosity function is with the functional form A. This is the LF of the spatially homogeneous 
sample. 
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Fig. 17. — The input and estimated LFs of the mock 2dF catalog. The cosmological parameters 
are Oq = 1.0, Aq = 0, F = 1, and as = 0.55. 

Fig. 18. — The same as Figure 17, except that the cosmological parameters are JIq = 0-3, Aq = 0.7, 
r = 0.25, and as = 1.13. 
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Fig. 19. — The same as Figure 17, except that the cosmological parameters are Qq = 0.3, Aq = 0, 
F = 0.25, and ag = 1.13. 

Fig. 20. — The /-band LF derived from the HDF photometric redshift catalog prepared by 
Fernandez-Soto et al. (1999). The dotted line represents the local /-band LF obtained by Metcalfe 
et al. (1998). 

Fig. 21. — The same as Figure 20, except that this is the S-band LF. The dotted lines depict the 
B-band LFs derived by Sawicki et al. (1997). 

Fig. 22.— The B-band LF of the HDF at the redshift range < z < 0.5. The dotted line depicts 
the Sawicki et al. (1997) LF, dot-dashed line represents the Metcalfe et al. (1998) 5-band LF, and 
long dashed line is the Autofib LF at z < 0.1. Our LF and that of Sawicki et al. (1997) agree 
with higher-normalization LF reported by Autofib Survey, but are significantly higher than that of 
Metcalfe et al. (1998). 

Fig. 23. — The evolution of the luminosity density, pl. Upper panel: the B-band luminosity 

density. Lower panel: the /-band luminosity density. Open squares are Ph{B) derived from CFRS 
(Lilly ct al. 1995), open circles, Ph{B) from Autofib (Ellis ct al. 1996), open triangle represents 
the value from Stromlo-APM (Loveday et al. 1992) and open diamond, ESP value (Zucca et al. 
1997). Crosses are the estimates of Sawicki et al. (1997). In this paper we did not try to correct 
for the reddening effect of dust. 

Fig. 24. — The result of the polynomial fitting to the /C-correction of galaxies constructed from 
the data prepared by Kinney et al. (1996). Symbols represent the representative SEDs of galaxies 
(open squares: elliptical, open triangles: Sbc, diamonds: Scd, and crosses: Irr). Panels (a) - (f) 
correspond to the fitting order 1st - 6th, respectively. 
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Table 1: The AIC of the polynomial fitting model for iiT-correction 



Type 




2ud 


Polynomial fitting order 

3rd ilh 


51 h 


Clh 


Elliptical 


96.4724 


41.6174 


42.9716 


19.0966 


-6.93137 


2.50965 


Sbc 


63.4433 


-9.27677 


-7.30190 


-26.7273 


-45.7942 


-32.1221 


Scd 


5.04120 


4.66447 


-0.536911 


-12.5160 


-11.0097 


-24.8293 


Irr 


-7.86433 


-31.3419 


-31.0330 


-49.8901 


-65.2359 


-57.3888 



